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ABSTRACT 


The Finite Element Analysis Program (FEAP) was expanded 
to solve linear and nonlinear, two and three’ dimensional 
heat conduction problems. The usual types’ of. boundary 
conditions, including radiation, may be specified. A wide 
range of two- and three-time level schemes for the solution 
of time dependent problems is available in the program and a 
discussion of those most commonly used is presented. The 
algorithms for the solution of typical practical problems are 
described and several numerical examples are presented. The 
results are compared with the available analytical 
solutions. A listing of this expanded version of FEAP and 


the corresponding user's instructions are provided. 
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1. INTRODUCTION 


This thesis describes an application of the Finite 
Element method to Heat Transfer analysis through the use of 
a computer program. This program was designed to solve 
linear and nonlinear, steady and unsteady, two and three 
dimensional heat conduction problems’ involving temperature 
dependent thermophysical properties and complicated 
radiation/convection boundary conditions. 

The Finite Element Analysis Program (FEAP), programmed 
by Professor R. L. Taylor in the Department of Civil 
Engineering of the University of California, Berkeley, 
served as. the point of departure for the present code. It 
has now been’ expanded with two additional modules for Heat 
Transfer analysis and time integration of first order 
equations, respectively, but the original characteristics of 
a research and educational tool in which the various modules 
can be changed or added to as desired, were maintained. 

Although all equations are retained in core, the program 
can handle realistic engineering problems with several 
hundred unknowns in most computer systems. 

While the algorithms to solve steady heat’ transfer 
problems are well studled and defined, more research has to 
be done on the solution of unsteady problems. The 
Zienkiewl cz two- and three-time level schemes were 


Integrated In the program In such a way that a wide range of 
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choices of time integration algorithms is available to the 
user. A limited study was performed in order’ to determine 
the characteristics of the most used members of those 
families of numerical schemes. 

The FEAP code is written in FORTRAN IV language and this 
version was constructed and tested on an IBM 360/67 computer 
system with 0S/67 release 18. Other systems and 


installations will require modifications. 
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11. STATEMENT OF OBJECTIVES 


A. EQUATION OF CONDUCTION. INITIAL AND BOUNDARY 
CONDITIONS 
The problem considered for solution fs thoroughly 
developed by Arpaci in Ref. 1 and is mathematicaly 
described in the region 2 by the equation 
pveie—agy {kK VI.) +O (1) 


subjected to boundary conditions 


T=T on Ty (2) 
Ww 
and 
CokeVE en, slg + quer og. 0 on [Ts (3) 
and initial condition 
lim T = T (4) 
t+0 0 


Vis the gradient operator, 1, and T2 are mutually 
exclusive parts of the boundary of the region 2 , T is the 
temperature and t is the time. The thermal eapaci ty oc, the 
thermal conductivity K and the rate of internal heat 
generation Q are thermophysical properties dependent on 
temperature, 

In equation (3) HF is an outward unit vector normal to 


the boundary surface whlle q, q represent respectively 


Cc? qy 
the Imposed heat flux and the rates of heat flow per unit 


area due to convection and radiation defined as 


geeghe CAT eT) | (5) 


and 


: a Ss gl ne . 2s _ i 
' t ia - Fr 7 7 
a a if rt soe 7 0 


Tie Wares 


a | So) Pore oe ta ae 


yisauavon? 2! otsiitos wor Tepe rr ong 
vinnizemerzgem ci tne f taf af  Ipegth ¥ sie 

noliaupe edy ed b no tye nz “7? 
(x3 9° 0T) Ve ae 


angis lbaao ereboued a 


wa i. 
ox p* no ; eM : 


= 
(ft) 7 mo 0p fp hae 

nod Vonage 4 
‘a3 at » T. 
vilewdum e676 ci Sree :1 .%oveveqo Seelbarg t 
wits et 7 2 nolnet of3 to yeebnudd of3 Yo emt 


at’ 99 yllasaso tamiedt? off amis etd el 3 Baas 


‘oad tnnwetol Yo e363 off bam R eeivigs 
no Snebneqeb eal y4eqora ac tong 9 a8 Hes 


va : it 7 4 4 

; 7 
03 (amen Yolzav Dinu baswsle ne 2) 7 (2)... no 
yievitssaze? 2ne2s 1907 a? ba neal 
tinu seq wet? ised 19 stoi 2 


7Bu 


4 


4 
qd. \* Eo #@ Teoe Toe 


Vga h,. ( T - Ta ) (6) 

In (5) he is the temperature dependent convective heat 
transfer coefficient. In (6) the parameter h. is defined by 
the expression 


h.=Fo (T¢ + 7 


< air Cet + Tie ) , where F is 


the radiant exchange factor ando the Stefan-Boltzman 


constant. T and Tae are the equilibrium temperatures for 


ac 


which, respectively, no convection and radiation occurs. 


B. NUMERICAL SOLUTION. DISCRETIZATION BY THE GALERKIN 
METHOD. MATRIX FORMULATION 
The spacewlse discretization of equation (1) using 
cartesian coordinates can be acomplished by Galerkin's 
principle as shown by Zienkiewicz in Ref. 2 and Lew in Ref. 
Bis 
Let the unknown function T be approximated, throughout 
the solution domain at any time t, by the relationship 
Z = EAN Coy,2) T(t) = <N> CEE (7) 
where N. are the usual shape functions defined piecewise 
element by element, T; being the nodal parameters. 
The result is 
[K] {T} + [c] {T} + {F} = {0} (8) 
where [K] and [C] are symmetric matrices defined, on the 


element level, as 


T T 
em éN. _6N SN. _6N 
[K]° = Ie (ax? — k,, + Sas Sones Ls + 
6N ON 
<ar <3? k,) dQ + 
Jie ( he + hy ) <N> <N> dI (9) 
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The load vector F is 
(Fy? = -[ ous @ davt 


vent 
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ole yar (11) 
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where the surface integrals are performed only over’ the 
surfaces where the prescribed boundary condition applies. 

It must be noted that the set of equations (8) is 
nonlinear since the matrices [K] , [C] and vector {Flare 
dependent on T. 

The system of equations (8) may now be solved by any 
automated numerical technique. Code FEAP was written for 


this purpose and Is described in the rest of this thesis. 
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Ill. EMNITE ELEMENT ANALYSIS PROGRAM 


The Finite Element Analysis Program (FEAP) was 
programmed and published by Taylor in Ref. 4 where it is 
well explained. The reader should consult this reference for 
details about it. 

This thesis is an expansion of the original program 
towards the solution of Heat Conduction problems and as such 
only the main features of the basic program will be referred 
in this section. A more detailed explanation is reserved for 
the new subroutines added to the program. The complete 
user's Instructions for FEAP are given in Appendix A and 
they complement this section. 

A. MAIN PROGRAM 

FEAP can be separated into two basic parts: 

a) Data input module and preprocessor 

b) Solution and output modules 

The data input module must transmit sufficient 
information to the other modules so that each problem can be 
solved. All the input data [fs stored in a single array 
(integer array M) which is partitioned to store all the data 
arrays, as well as the global arrays, e.g., stiffness, mass, 
load, etc. 

The total capacity of the program Is controlled by the 


dimension of the array M tn the’ blank common of the main 
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program and the corresponding value of the variable MAX. 

The partition of M Is performed in the control routine 
(subroutine PCONTR) and the data used ffor it is supplied by 
the user in’ the title and control information cards (first 
two cards for every run of the program). 

Then the nodal coordinates, element connections, 
material properties, nodal loading and boundary conditions 
may be input using the macro control statements. They are 
interpreted by subroutine PMESH and each control is a 
function independent of the others. 

To clarify what was said till this point, the discussion 
of a simple example may be useful. Consider the mesh 
presented in Figure l and assume that quadratic rectangular 
elements as the one defined in Figure 2 are used. This 
example is completly solved in Appendix B and will be used 
throughout this text. 

The sector of a circular ring was divided in six 
elements and the 33 nodes were numbered. The order of 
numbering is not crucial but in order to improve the profile 
of non-zero coefficients the following general rule should 
be used. 

The numbering should be such as to minimize the nodal 
difference for each element (maximum node nummber’ minus 
minimum node number). 

The user can now proceed to the preparation of data for 
the program. The first step consists of specifying the 


problem title and the control Information, The latter Is 


‘Figures are grouped at the end pp 47-79 
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clearly 33 nodes, six elements, one material set, two 
spatial dimensions and eight nodes per element. If this is a 
Heat Transfer problem we have one- unknown per node 
(temperature). 

The program expects now the data cards for mesh 
description. Any analysis require at least 

a) coord|nate data which follows the macro command COOR 

b) connectivity table which follows the macro command 
ELEM 

c) material data which follows the macro command MATE 

In addition, most analysis will require specification of 
nodal boundary restraint conditions, macro BOUN, and the 
corresponding nodal force or displacement values, macro 
FORC. The term displacement refers to the value of the 
unknown for each specific problem. In a Heat’ Transfer 
problem the unknown Is temperature while in’ an Elasticity 
analysis it is a geometrical displacement. 

The TEMP macro command (temperature data) shall not be 
used In a Heat Transfer analysis to specify nodal 
temperature. It may be used to specify auxilary nodal 
quantities. 

If in our problem we assume the line defined by nodes 
31, 32 and 33 at the temperature of 20, the line defined by 
nodes 1, 2 and 3 at 820 and Insulated elsewhere, the macro 
control statements and correspondent data cards will be 


COOR 
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31 a) 
4 5.116666667 


29 - 283333333 
2 Sate L 3. 
32 3 36 
3 Sire 6 
33 3 6 


5 5.116666667 6. 
30 32835553555. 6. 
(blank card) 
POLA 
1 33 1 
(blank card) 
ELEM 
1 1 iL 6 8 3 
(blank card) 
BOUN 


33 -1 
(blank card) 
FORC 


31 i 20. 
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(blank card) 
MATE 


10. 10. 8000. 250. 3 
END 

In spite of the fact that the MATE command must be 
included in this part of the data preparation, its 
discussion is left to the section where the element modules 
are treated, 

After completing the mesh data input, we are’ ready to 
Initiate the problem solution. 

FEAP has modules for variable algorithm capabilities and 
which, If necessary, can be modified or expanded. The basic 
aspect of the variable algorithm program is a macro 
Instruction language which can be used to construct modules 
for specific algorithms as needed. This language is 
Interpreted by subroutine PMACR and the complete list of 
macro instructions available in this version of FEAP is 
given in Appendix A.3. 

Here, as an example, we will discuss in detail the 
algorithm for the solution of linear steady state problems 
using the original explanation given in Ref. 4. 

To use the macro programming commands, the user. only 
needs to learn the mnemonics of the language. If one wishes 
to form the global stiffness matrix the program instruction 
TANG Is used (TANG Is the mnemonic for a symmetric tangent 


stiffness matrix and for nonlinear elements would form and 
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assemble into the global stiffness the element’ tangent 
stiffness computed about the current displacement state; for 
linear elements this is just the linear stiffness matrix). 
For a problem with an unsymmetric tangent stiffness the 
macro command UTAN Its used. 

If one wishes to form the right hand side of the 
equations modified for specified displacements one uses the 
program instruction FORM, The resulting equations are 
solved using the instruction SOLV. 

Printed output can be obtained using the macro command 
DISP, 

The above instructions are sufficient to solve linear 
steady state problems, that is, the macro [Instructions 
TANG ( or UTAN ) 

FORM 

SOLV 

DISP 

are precisely the required Instructions to solve any linear 
steady state problem. 

If the same block of instructions is repeated twice or 
more, considerable effort is wasted in preparing the macro 
instruction data, To rectify this, looping commands are 
introduced as the Instruction pair 


LOOP n 


NEXT 


which tndicates that looping over all instructions between 
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LOOP and NEXT will occur n times. 

Many other classes of problems can be solved using the 
macro instruction list given in Appendix A. With the above 
short explanation and the discussion presented in section IV 
for Heat Transfer problems, the reader will be able to 
construct his own algorithms using the macro.- programming 


language. 


B. FINITE ELEMENT SOLUTION MODULES 

Most of FEAP [fs common to a wide range of problems, but 
some of the computations are linked to the’ type of problem 
to be solved. It is clear that if one wishes to’ solve an 
Elasticity problem, the calculation of the stiffness matrix 
will be different from that used for a Heat’ Transfer 
problem, It is also clear that differences also arise if one 
is using triangular elements instead of quadratic 
isoparametric elements as a way of describing the continuum. 

The program was designed such that all computations 
associated with any type of element are contained in an 
element subroutine called ELMTnn where nn is between 01 and 
05 in this version of FEAP. Each element type to be used is 
specified as part of the material property data, following 
macro control MATE. 

Since the objective of this work is the solution of Heat 
Transfer problems, the element modules discussed here are 
the two and three dimensional heat tranfer modules named 


ELMT02 and ELMTO3, 


18 


at ener al (hu i : 

| v4 ia Be aaa " 
ert artheu bevior od na) ‘bmatdowa * 
svods one utgtWl A xi bnaagaent nevis 3 
Vi cotgsoe nt besnerer eplaewel> odd fap” 
of sfds ed fTlw sebsen od .emetdotg 
En lemetgte ovoem od ‘grlen emt stron t® nwo €h 

>) eae 


ae os 


ez1v00M woITUJOe THEME: 
jud .ametderq Yo exneay sbhiw 6 0% nonmos «! 9433 to 


na’ 


io%d Yo gays sAJ oF bavntl ave anc tyaguqnes — 
ne OVige oF eodziw eno 7) sadd seeld ef Bf 
zixjem seentttiae ef2 YO neolzaluelag saz snote 
yetenatT jee #® Oo? bewd get? “manl treceT 
: 

no %1 aziz cel e209 neve Tt th fend yeals oefe 204 

ot jeubéeup Yo bessenl ashemate valognelye 


smuuatano9 oft anidtsoeeh Yo wer & 46 adnemets So 
ar 


enoltatudmon fife ted? oue bengtesh ow megs 
os of fentletnos st6 Inanels Foe Saya vhs a3 For 
bie 20 aeewted el fn eNniw wet Ad 3 belias onls 7161 
Neau sd of sava snemete toed W835 Yo aovetay ‘ 
gitwotfo? (ssab vdseeetq fsalvetan ans Yo 218q 
3TAM to 
$604 Yo weltuloe oft at anew Bias to ov isaetdo, ) 
ew ooo feeouselb goluborm Items ts ony 


bemat eeluben vetnbss 
* 


ee 


1. Iwo Dimensional Heat Transfer Element 

Subroutine ELMT02 is accessed using the number two 
in column ten of each material number card. 

According to the value of the variable ISW, defined 
in subroutine PMACR, a specific function ts performed. If 
1SW=1 it reads and prints the material property data and 
line boundary conditions. The stiffnes and mass matrices are 
computed if ISW=3 and ISW=5 respectively. The load vector 
due to internal forces and boundary conditions is calculated 
when ISW=6. 

The material properties may be input as constants or 
as a table of temperature and property values. If any other 
way of defining the property values is required by a 
specific problem, subroutines PKX , PKY , PROC’ or PQ which 
calculate respectively the conductivity in x and y 
directions, the thermal capacity and the heat generation per 
unit volume, can be easily modified or substituted. 

The shape functions used by this heat tranfer module 
are calculated by subroutines SHAPE and SHAP2. These 
routines are capable of constructing shape functions for a 
three-node triangle, a four-node IJInear quadrilateral, an 
elght-node quadratic serendipity quadrilateral, a nine-node 
quadratic Lagrangian quadrilateral or any combination 
in-between. 

A four- to nine-node two dimensional element’ is 
shown in Figure 3. As shown by Bathe and Wilson in Ref. 5 


the temperature within the element is expressed at any time 
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in the local coordinate system s, t in terms of the nodal 
temperatures Tj by 


ks, tose 2 N.(s,t) T; 


where 

N, = (L-s)(1-t)/l - (Ng#Ng)/2 - No/4 
No = (L#s)(1-t)/4 - (No#Ng)/2 - No/t 
Nz = (1#s)(1+t)/6 - (N,#NQ)/2 - Ng/t 
Ny = (les) (1+t)/4 = (N74Ng)/2 = No/4 
No = (1-s7)(1-t)/2 

Ne = (1+s)(1-t7)/2 

No (1-57) (1+t)/2 

Ng = (1-s)(1-t7)/2 

Ny = (1-s2)(1-t2) 


If any of the nodes from five to nine are omitted 
the corresponding value of N is zero. 

The connectivity table must follow the numbering 
convention shown in Figure 3, otherwise wrong results will 
be obtained. 

The maximum number of nodes’ per element will depend 
on the type of elements to be used and may be from three to 
nine. The eight-node serendipity shape functions occur if 
the ninth node number is omitted. The four-node 
quadrilateral shape functions are computed if only the first 
four nodal connections are non zero and the’ three-node 
triangle if the first three nodal connections are non zero. 
Finally, if a mid-side nodal connection is omitted the edge 


Is linear. 
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The line boundary conditions, closely related to the 
type of element used, are defined in the cards following the 
macro MATE. The line subjected to a specified boundary 
condition is identified using the local coordinates. A line 
is numbered 1 or 2 if it is perpendicular to the axis s or 
t, respectively. Those numbers’ are positive or negative 
according to which direction of the axis it is 
perpendicular. Thus, as an example, the line defined by the 
nodes 2, 6 an 3 Is numbered 1 while the line defined by 
nodes 1, 5 and 2 is numbered -2. 

The boundary conditions are treated in subroutine 
BCOND2 and four types are considered: 

a) specified flux 

b) convection with constant heat tranfer coefficient 

c) convection with temperature dependent heat 
transfer coefficient 

d) radiation 

If more complex boundary conditions are’ required 
subroutine BCOND2 may be easily modified. 

The numerical integration necessary to evaluate the 
matrices and vectors in equation (6) are performed using the 
Gauss quadrature formula; the abscissas and weight 
coefficients are calculated in subroutine PGAUSS. The user 
may choose the number of points’ per direction used in the 
integration from one to six, but the default value is four 
points , 


To illustrate the above, we assume that our sector 
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used as an example is part of the cross section of a_ hollow 
cylinder made of a material with conductivity in x andy 
directions of 10, specifle heat 250 and density 8000. If the 
outside surface is subjected to convection to an ambient 
atmosphere at 20 degrees with a constant heat’ tranfer 


coefficient of 200, the cards following the macro card MATE 


are : 

MATE 
L 2 

10. 10. 8000. 250. 3 
6 2 1 200. 20. 


It must be noted that a plane analysis is specified 
and the integration will be performed using three points per 
direction. The reader must also note that this’ version of 
the problem is slightly different from the one’ given when 
the mesh data preparation was discussed. In the’ previous 
problem the line now subjected to a convection boundary 
condition was at a specified constant temperature of 20. If 
this second problem is to be solved the cards corresponding 
to nodes 31-33 in the BOUN and FORC sets must be omitted. 

Following is a_ list of subroutines included in the 
two dimensional heat tranfer module: 

ELMTO2 : main routine; forms the necessary arrays 
for the solution, 

SHAPE : calculates the shape functions for. the 
triangle and the l!near quadrilateral elements. 


SHAP2 : adds the quadratic terms and the center node 
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to the shape functions. 

PGAUSS : determine the abscissae and weight 
coefficients for the numerical integration. 

BCOND2 : adds to the arrays calculated in ELMT02 the 
contribuition from the specified boundary conditions. 

PKX : determine the temperature dependent 
conductivity in the x direction. 

PKY : determine the temperature dependent 
conductivity In the y direction. — 

PROC : determine the temperature dependent thermal 
capacity. 

PQ : determine’ the temperature dependent heat 
generation per unlit volume. 

CONV : determine the temperature dependent heat 
transfer coefficient. 

TABLE : called by PKX, PKY, PQ, PROC and CONV; 
calculates’ the correspondent coefficient to a given 
temperature by linear interpolation between two consecutive 
entries in a table. 

JACBB2 : calculates the jacoblan determinant when 
an integration over a line is necessary. 
2. Ihree Dimensional Heat Transfer Element 

This module is called ELMT03 and is accessed using 
the number three in column ten of each material number card. 

This element is the generalization of ELMT02 for 
three dimensional space and little remains to be said about 


its use. 
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The material properties are read in a similar way as 
the two dimensional case, but the conductivity in the z 
direction must be added. This module’ uses the same 
subroutines as ELMT02 to read and calculate the temperature 
dependent properties. 

Subroutines SHAP3D and SHAP3 may construct 
interpolation functions for. any combination between a 
eight-node linear brick and a 21-node Lagrangian brick. 

An eight- to twenty one-node three dimensional solid 
element is shown in Figure 4&. The natural coordinates (r, s 
and t) of the eight corner nodes are (+1, +1, +1) , of the 
twelve mid-edge nodes are (0, +1, #1), (#1, 0, +1) and (+1, 
+1, 0) and of the center node are (0, 0, 0). 

The temperature within the element is defined in 
terms of the nodal temperatures T; at any time by [Ref. 5] 

T(r,s,t) = 2 Ni(r,s,t)T; 

If we define 

G(8,B;) = -5(1+8;8) » for B.= +1 
= 1-87 
and 
g. = G(r,r;)G(s,s;)G(t,t;) 
where rin Sq and t; are the natural coordinates of the 


element nodal points, the interpolation functions Ns are 


= (Ng +N 7/2 2 N,4/8 


Tea 


=g = (Ny +N #N,Q)/2 - No1/8 


2 10 
Nz = Bz 7 (Ny gtNyytNyo)/2 - Noy/8 
N 


pen ed megan tae var N29) / 28 N50/8 
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he = g; - N54/4 for j=9,...,20 
idle Sk 


If any of the nodes from nine to twenty one are 
omitted the corresponding value of g is zero. 

The numbering convention used in the definition of 
the shape functions and shown in Figure & must be followed 
by the user when’ the connectivity table Is bullt and the 
surface boundary conditions are specified. The surfaces are 
numbered +1, +2, +3 as they are perpendicular to the 
positive or negative directions of ne. axis r, s, t, 
respectively. 

As the numerical integration routine is the same as 
for ELMT02, the user may choose from one to’ six points per 
direction, the default value being four points. 

Following is a list of subroutines included in the 
three dimensional heat transfer module: | 

ELMTO03 : main routine; forms the matrices necessary 
for the solution. 

SHAP3D : calculates the shape functions for’ the 
elght node linear brick. 

SHAP3 : adds the quadratic terms to the _ shape 
functions and the center node for the Lagrangian brick, 


BCOND3 : adds to the matrices the contribuitlon from 
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the specifled boundary conditions. 
PKZ : determine the conductivity in the z direction 
JACBB3 : calculates the jacobian determinant when an 


Integration over a surface is necessary. 


C. TIME INTEGRATION MODULE 

The first order ordinary differential equation solver is 
accessed when the macro command ODE1 is utilized and uses 
the Zienkiewicz two- and three-time level schemes. The 
details of these algorithms are treated in section IV. 

This module was designed for the solution of first order 
equations but can, with Jlittle programming effort, be 
expanded in order to include higher order algorithms and 
solve higher order differential equations. 

In order to economize computer memory space at some cost 
of execution time, this module only uses the space 
previously reserved for the mesh construction. As the mesh 
data Is not needed during the execution of atime step 
integration, all the corresponding arrays are kept in file 9 
during that period. Once the next displacement vector. is 
calculated and the current displacement vector is' saved in 
file 10, all the Information in file 9 is retrieved and the 
mesh data restored. 

File 10%s used as a working space during the time 
integration and the current displacement vector is saved in 
It for possible future use. This vector is retrieved if the 
algorithm used in the next time step integration is a three 


point scheme, 
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Besides the time step size change using the appropriate 
macro instructions (macro TOL), an optional automatic time 
step adjustment was incorporated in the program. The norm of 
the difference between the displacements vectors at two 
consecutive times is computed at each step. If the norm is 


less than a_ predetermined value AT the time step size 


max ’ 
is doubled before going to the next step, whereas if the 


norm is greater than a_ prespecified value AT AG the time 


in’ 
step size is halved and the calculation for that time step 
is repeated until the norm is acceptable. The magnitudes of 
the maximum and minimum values of the norm are’ problem 
dependent and must’ be choosen by the user. If they are not 
specified no time step adjustment will be performed. The 
recalculation of the displacement vector is allways done 
using the two. point scheme, even when a_ three point scheme 
was utilized for the first calculation. 

The macro command ODE1 used to access’ this module 
must be followed by a second macro command in order’ to 
determine the function to be performed. The secondary macro 
instructions are INIT, LINE or QUAD. 

The couple ODE1 INIT reads the integration constants 
theta , beta and gamma, the parameters for the automatic 
time step adjustment and the [Initial displacement vector. 
This data must follow the macro program (see user's 
Instructions). No tIlme integration is performed by this 


Instruction, 


The couple ODE1 LINE performs the two point scheme 
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and the current displacement vector is substituted by the 
new calculated displacement vector. 

The couple ODE1 QUAD has a function similar to ODE1 
LINE but uses the three point scheme. 

Once the execution of an ODE1 macro command is 
complete, the mesh data is restored and the displacement at 
time t replaced by the displacement at time t+ At, when the 
secondary macro command is LINE or QUAD. If INIT is used, 
the displacement vector becomes the given initial 
displacement vector. 

The subroutines included in this module are: 

PODE1 : control routine; reserves working space and 
calls the appropriate subroutines. 

INIT : reads the constants theta, beta and gamma and 
the maximum and minimum values for the norm used in the time 
step adjustment. 

PLINE : performs the two point integration scheme. 

PQUAD : performs the three point integration scheme. 


TERM : performs the automatic time step adjustment. 
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IV. APPLICATION OF FEAP TO HEAT TRANFER PROBTEMS 


The program FEAP requires from the user the knowledge of 
the algorithm to be used in the solution of the problem to 
be treated. 

In this section the algorithms used in’ the solution of 


heat conduction problems are discussed. 


A. STEADY STATE PROBLEMS 
1. Algorithms for Linear Problems with Macro Program 


In this case equation (8) becomes 


[Kk] {7} + {F} = {0} (12) 
The algorithm used to solve (12) is described by 

[Kk] Avi se={F}°=(K] 11°F (13) 

(eit. iv} (14) 


where {T°} are the specified nodal temperatures. 

The macro instruction TANG builds the left hand side 
of (13) while FORM forms the vector in the right hand side. 
The instruction SOLV solves the system (13) and performs the 
step defined by (14). 

Consequently the simplest macro program used _ to 
solve thls problem Is 
TANG | 
FORM 
SOLV 
DISP 
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2. Algorithms for Nonlinear Problems with Macro 
Program 
Equation (12) still applies to this case but now [K] 
and {F} may be dependent on’ the nodal temperatures {T}. The 
well known Newton-Raphson iteration may be used and _ the 
algorithm [ts described as 
[kGlt Wetveia= -CFET*)} - [keT*}] {T*} = (Ry 025) 
(Te eat =A THe eat) (16) 
where {T1} are the nodal temperatures at the ith iteration. 
The algorithm defined by (15) and (16) is the one used for 
the linear problem repeated several times. Then’ the macro 
program may be 
LOOP n 
TANG 
FORM 
SOLV 
DISP 
NEXT 
DISP 
where n are the user's guess of the number of iterations 


necessary to obtain equilfbrilum. However, the program has an 


internal check on the value of the vector norm [Ir*|| 7 
where 
a eal/ 2 
HIRI] = @ [REI 
k=1 
Whenever 
[|Ri|[< ToL x max [IRI || Cite...) 


J a 
where TOL Is a predefined tolerance ( 10° ts the default 
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value), the iteration ceases and a skip to the macro command 
immediately following the first NEXT occurs. Usually one 
begins with zero as the initial guess of the displacement 


vector; however, any other vector may be used. 


B. TIME DEPENDENT PROBLEMS 
1. Two and Three Point Recurrence Schemes for First 

Order Equation 

The set of differential equations 

[K]{T} + [C]{T} + {F} = {0} (8) 
may be solved using one of the many recursive schemes 
described in the literature. From them, the two. and three 
time level schemes presented by Zienklewicz [Ref.4] offer a 
wide range of choices for the solution of linear and 
nonlinear problems, 

The recurrence relation for the two-time level 
scheme can be written as 

( aelC] + OLK]){T,,,} + (-gelC] + (1-0) [K]){T,} + 

{F} = {0} ; O7< 0 <1 (17) 
where the subscripts mn and n+l denote evaluation at time t 
and t+at. The vector {F} Is defined as 


{F} = {F..,} + (1-6) {F} (18) 


n+l 

The choice of the parameter © defines the particular 
scheme to be used and the reader will recognize a well known 
series of finite difference formulas with a modification of 
using a weighted loading term {F}. 


Consider the system of decoupled equations in terms 


of the modal participation varlables Tj 
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C5 T; + k; T; + f. = 0 (19) 


for free response, i.e., f;=0, expression (17) becomes 
1 1 = 
Grete KO) AU ee miGae cy) + ki 9)) (T,),,70. 20) 


Substituting the relation 


(Tipe 


in (20), the resulting characteristic equation of the 


=H ae (21) 


recurrent scheme solved for i) gives 
A = [1-k; (1-0) At/c;]/[1+k,0 At/c,] 
The scheme will be unconditionally stable if 
[ral <j 
for any value of P; where 
Pia (k;/c;) At (22) 
This condition is satisfied for 
Q.2 1/2 
On the other hand if 
O26. 160'< 2 
stability is conditional requlring 
Py < 2-26) 

Figure 5 shows how A varies’ with P; for the schemes 
discussed later and how {ft compares with the exact value of 
d = exp(-p;) 

The schemes considered are 

a) Crank-Nicolson with 0=1/2 

b) ©=2/3 proposed by Zienkiewicz [Ref. 4] 
c) 0=3/4 proposed by the author 

d) ©=0.878 proposed by Liniger In Ref. 6, 


The Zienkiewicz three-time level scheme is defined 


32 


(er) ; er 


 geméoed ctr? sa ot ar * = 


(os) 0=,f,7) (CO pede pa ad 
nulbetes of39 nite 

(15) agtyTA- 
efit .%o nplieupe Vgelvesoetaiio~ galt luest 
roviss A yo" bevios 


ah Oy IML orate f) ae 


¢ 
+] el@are vi swo)o Minoo ed Tih iw 


os arte 49 ‘0 
ceo) 16 €,a\34) “3 


— 
ee 


9 bealvetdce ol noi? lhnog 


a (a r 


if ’ 


Y) boat senze 


gn '{*iupss (anol? Tonge 
- 


a 


Fos DP) Xi 2 

sometioe ony 20? mo ad'4 zelnay 4 wort ewole @ 
Yo auilev JDGKS. end, Hai” schon 3) wert bas 
(,a-)ae * 

pve wotab ien0d sseaae 

Silee fzlw ooe lenis 

[a teh) eaiwetanatl wd & 

ae ” edt 


| 


as 
(y[C] +8 At[K]){T,,,} + ((1-2y)[C] + (1/2-28+y) At [K]) {T,}+ 
(-(1-y) [C] + (1/2+8-y)At[K]){T, 4} HONCLE T2240 (23) 
where the subscripts n, ntl and n-1 denote evaluation at 
time t, t+At and t-At respectively. 

The vector { F} is defined as 
Ama RE ee +L /2-2ery te} eC /2te-yItF! 4) (24) 
and the parameters y and g define the particular scheme to 
be used. 

Considering again the system of decoupled equations 
(19) and using the relations (21) and (22) one’ finds that 
the characteristic equation for the three-time level scheme 
is 


(y+ep;) A” 


+ [(1-2y) + (1/2-28+y)p,J0 + 

[-(1-y) + (1/2+8-y)p,] = 0 (25) 
Writing 
g = [1+(1/2+y)p,]/Ly+6P,] 


h = [-1+(1/2-y)p,]/[y+8p;] 


the roots of (25) may be written as 

Agig = Ree) /2e [(2-g)7-4(1-ny]4/2/2 

These roots will be complex if the quantity under 
the square root is neeatioe. Then the modulus of i is 

lalie engi? 

The scheme is stable if |A|<l for any p,and Wood in 
Ref. 7 shows that this condition is satisfied if 


y>1/2 and B>y/2 


33 


- coy aa is ' —— 
yr : mi if ‘7 “ By Ah 
ir 4 , _ a 


* oper i hee: a 
at eee ses aivs-49) a _ 
ces)! OY ae sips its rm | fs 
so noliauleve sidneb tes ‘bas fen ar ‘i708 2. 


xtovi dosqaet 3 hed br 
“a8 bon! ete: at 


ceo) 2° Cy gh) Cer mea hed { aTiteeasteysy. we 
oi smatise veluolores of2 an tho a bre ¥ 23 


enolreune teatauesed Yo mateve ofa alesse aninee ? nM 
refiz ebal? sno (St) bre CLD) anolastew ats ante 


gnatioe tevel emit-osrd? eff ww golisuns side 


hE aye S\0), > (yS-L)) i 
vesy 0 = [,abyserSNEh w tye 8)-3 . P 
; a a 
t ahew INT pals S\ 1) +S ae 
ati “\t) a 
: 
ae reap toh alincao east 
eying. ye" (gS a¥te-t) *« 
sébnu valinevp ef7 11 x9 Teta “pa i j atees 
2e| 4 Yo eutubom edd mubates. 


_ 


; m es aymiA 


_#} bool bane ves OF sii ee 
& boiNal: i 


Equation (25) has two roots: ho , the principal root 
which is the aproximation to the exact 
d= exp(-p,) 
and the spurlous root os For relative stability one must 
have 


2! 


Also a negative root or complex roots can produce 


heen 


oscillation. 


1 and ho and the 


modulus |\| for the:complex roots are plotted against Ps 


In Figure 6a-f the real roots i 


for the schemes: 

a) 8=1/3, y=1/2 proposed by Lees in Ref, 8 

b) g=3/4, y=1 proposed by Hogge in Ref.9 

c) 620.646, y=1.2184 proposed by Wood [Ref. 7] 

d) 6=4/5, y=3/2 proposed by Zienklewicz [Ref. 4] 

e) B=9/10, y=3/2 proposed by the author 

f) Fully implicit algorithm, §=1, y=3/2 [Ref. 7]. 

From Figure 6a one may predict a very” strong 
oscillatory behaviour for the Lees algorithm a). 

In order to study the various schemes, the physical 
problem represented by the nondimensionalized linear heat 


conduction equation 


éT 677 
ot 6x* 


was solved in a bar of length 4& and width 1 subjected to the 


boundary conditions 


T=1 at x=0 
ST_-9060 at si=24 
6x 
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and initial condition 

; T=0 at t=0 
i.e., consider a step change in the surface temperature, 
This problem was also studied by Wood and Lewis in Ref. 10. 

The problem was discretized in space using ten equal 
length two dimensional quadratic elements through the length 
of the bar. This type of elements was chosen because it was 
reported by Wood and Lewis as’ giving the worst numerical 
results. 

In order to avoid the discontinulty in the loading 
term caused by the step change of the surface temperature at 
the Inittal time, the problem was first solved starting from 
time t=At, where At Is the time step size. The value of the 
temperature distribution at t=At iis provided by’ the exact 
analytical solution. This procedure also provides’ the 
necessary two starting vectors for the three-time level 
schemes. 

The temperature of the nodes at x=l is used as 
reference value to compare the various schemes. 

In Figure 7a-d the results obtained from the 
various two-time level schemes with At=2 for the temperature 
at x=l are compared with the analytical solution. For this 
Ideal starting conditions the Crank-Nicolson, 6=1/2, proves 
to be the most accurate algorithm. When the step size was 
Increased to 10, al] the schemes’ performed well as shown in 
Figure 8a-d. 


The corresponding results for the three-time level 
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schemes are presented in Figure 9a-f and Figure 10a-f for 
At=2 and At=10, respectively. 

Strong oscilations were observed with the Lees 
algorithm a). 

For At=2 the schemes c), d) and e) performed well. 
For the larger step size of 10 the safest schemes are e) and 
fi). 

The step change In the loading term has_ been 
reported (Ref. 10] as producing relative instability when 
some of the schemes are used, particularly with the 
Crank=-Nicolson algorithm. The noise may be  supressed 
reducing the magnitude of the time step but this is 
impraticable In most problems. Other techniques may be used 
to reduce the amplitude of the oscillations as_ those 
discussed by Wood and Lewis [Ref. 10] and Gresho and Lee in 
Ref, 11. 

The step change in the loading term in the problem 
in study [s due to the variation of the surface temperature 
at the initial time t=0. That discontinuity In the force 
term Is illustrated in Figure lla. When’ the time dimension 
is discretized It has been common practice to calculate the 
numerical solution at the nodes t=0, t=At, etc. When a 
three-time level scheme is used the node at t=- At is 
considered assuming all conditions steady before t=0. The 
force values used by this method are indicated by the 
circles In Figure lla, This procedure transfers to. the 


numerical solutlon the uncertainty In the value of the force 
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term at time t=0 introduced by the mathematical 
discontinuity assumed in the analytical solution. This 
problem is avoided if one uses the procedure illustrated in 
Figure 11b. The nodes considered for the discretization are 
those at t=-At/2, t=At/2, t=3At/2, etc. The discontinuity at 
t=0 is smoothed by the interpolation of the loading term 
inherent to the partIcular scheme being used. The values of 
the force term for every node are well defined and no 
uncertainty is associated with any of them. For’ the 
three-time level schemes the node at t=-3At/2 is used 
assuming all conditions steady before t=0. 

The problem studied was solved using both starting 
procedures, The results obtained with the first procedure, 
start at t=0, using the two- and three-time level schemes 
are shown in Figures 12a-d and 1l3a-f, respectively. The 
second procedure, start at t=- At/2, gives the results shown 
tn Figures lha-d and 15a-f for the two- and three-time level 
schemes, respectively. 

One may conclude that, In this simple problem, when 
the procedure Illustrated in Figure 11b is followed, the 
step change in the loading term do not deteriorate the 
accuracy of the numerical results obtalned with the schemes 
tested, with the exception of the Crank-Nicolson algorithm. 
This scheme shows an oscillatory behaviour not observed when 
the step change In the force term Is not present. 

The starting procedure of Figure 11b also proves to 


be an efficient way of starting the time Integration with 
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the three-time level schemes. 


2. Algorithms for Time Dependent Problems with Macro 
Program 


A simplified form of the two-time level scheme is 
avallable in FEAP and is described by 
(qeC,] + O(K,]){T 44} + (-elC,] + (140) (K,]){T,}+ 
KE{F} = {0} (26) 
where the subscripts n and n+l denote evaluation at time t 
and t+At, respectively. 

The three-time level scheme programmed in FEAP is 
described by the expression 
(vIC,] +8 at{K,1){7,,,} + CCl-2y)(C,] + (1/2-28+y) At [K,]) (7, } 
+ (-(1-y)(C,] + (1/2+8-y) 4tK,]){T,_1} 

+ te. } = {0} (27) 
where the subscripts n, n+l and n-1 denote evaluation at 
times t,At+ t and t-At respectively. 

It must be noted that jin (26) and (27) the vector 
{F tis an aproximation of the interpolated force vector {F}. 

If the force vector is strongly dependent on time, 
this aproximation may produce wrong results. 

In problems where the stiffness matrix is constant 
and the force term results from the specified boundary 
displacements or where the specified boundary forces are 
only tIlme dependent, the interpolation of the force term may 
be easily acomplished since the force/displacement boudary 
values can be changed at any time using the macro commands 


MESH or PROP. This procedure was found partcularly useful 
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when a_ step change in boundary displacements or forces Is 
applied at the initial time, 

The first order ordinary differential equation 
solver module is acessed using the macro command ODE1 
followed by one of the following macro commands: 

INIT to specify the initial displacement vector 

and the integration constants 

LINE to perform the two-time level algorithm 

QUAD to perform the three-time level algorithm 

Stnce the matrices [K] and [C] and the vector {F} 
must be formed before the time Integration, the macro 
commands TANG, CMAS or LMAS' and FORM are closely associated 
with the use of ODE1. It should be noted that the macro FORM 
forms the vector 

{R(T} } = aE 3 - [Kd 
which Is allways dependent on the current displacement. 
Therefore the macro command FORM must be used every time 
step and preceed ODE1. 

To solve a_ fully nonlinear problem the’ following 


macro program may be used: 


DT At 
ODE1 INIT 
LOOP n 
TANG 


CMAS ( or LMAS ) 
FORM 
ODE] LINE ( or QUAD ) 
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TIME 
DISP 
NEXT 

The macro command program must’ be followed by the 
data cards necessary to specify the integration constants 
6, y and 8 and the [nitial displacement vector. 

If the matrices [K] and/or [C] are not temperature 
or time dependent then the instructions TANG and/or CMAS or 
LMAS must be placed outside the loop. 

For the simplest case of a ]inear problem, the macro 
program may be 
DT At 
ODE1 INIT 
TANG 
CMAS ( or LMAS ) 

LOOP n 

FORM 

ODE1 LINE ( or QUAD ) 
TIME 

DISP 

NEXT 
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V. NUMERICAL EXAMPLES 


A. RADIAL TEMPERATURE DISTRIBUTION IN A HOLLOW 

CYLINDER 

Consider the hollow cylinder of infinite length already 
discussed in section I!1. Consider the case in which the 
inside wall is maintained at a constant temperature and heat 
is lost by convection through the outer wall. This problem 
was considered for solution by Lew [Ref. 3]. 


The geometrical and thermophysical characteristics are: 


Inside radius (r;_) 0.1m 
Outside radius (r.)4) 9.3 m 
Thermal conductivity (k) 10 W/m °C 
Outside ambient temperature 20 °C 
Heat transfer coefficient 200 w/m2°c 
Inside wall temperature (T5) 820 °C 


The finite element model of this problem was completly 
discussed before and iis shown in Figure 1. The complete 
computer output is given in Appendix B. It must’ be noted 
that the element arc width used is arbitrary since’ the 
problem is symmetric and no heat is conducted perpendicular 
to the radius. 

_ The comparision between’ the analytical and FEAP 
solutions is shown in Figure 15. In this case the values 


provided by the finite element model used are very accurate 


and no further mesh refinement is necessary. 
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B. NONLINEAR STEADY STATE HEAT CONDUCTION IN A SLAB 

WITH TEMPERATURE DEPENDENT CONDUCTIVITY 

A liquid is boiled by a flat electric heater plate of 
thickness 2L. The internal energy Q generated electrically 
may be assumed to be uniform. The boiling temperature of the 
liquid, corresponding to a specific pressure, is Twe 

In the finite element analysis for the temperature 
distribuition in the heater, ten equal length quadratic 
elements were used to represent the unit cross_ section 
through the half plate thickness L, since the problem is 
symmetric. The extremity of the resulting slab is assumed at 
the temperature Tg while the other is insulated. 

The thermal conductivity is assumed temperature 
dependent according to the expression 

k = a(1+bT) 

where a and b are constants. 

Assuming the following values for the parameters ina 


consistent system of units 


T = 0.0 
Ww 

a = 0.5 
L = 1.0 
Q= 1.0 


the problem was solved for b equal 0.0, 1.0 and 2.0. 
The numerical results are compared with the analytical 
solutions in Figure 16 and one may conclude that FEAP 


provides exact solutions in this problem. 
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C. TRANSIENT TEMPERATURE DISTRIBUITION IN A HOLLOW 

CYLINDER 

The same hollow cylinder of infinite length treated 
before is considered again but the case solved is the one 
with both outside and inside walls maintained at constant 
temperatures. Thermal diffusivity, a=k/ pc is specified as 
0.00005 m2/sec. Initially the cylinder is at an ambient 
temperature To of 20°C: suddenly the temperature of the 
Inside wall is raised to 820°C while the outside wall is 
maintained at 20°C. All other conditions are the same as for 
the steady state problem. 

In order to obtain the unsteady solution the number of 
elements in the finite element model was doubled, thus 
twelve radial quadratic elements were used. The element arc 
width chosen is the same 6° as before. 

The two time-level scheme with @9=2/3 was’ used for the 
time integration. In order to compare the results’ the 
dimensionless temperature T*, radius r* and Fourier number 
Fo were used. They are defined as 
elo ine ain 


T/Tin 


T* 


r* 
and 
Fo = at/ts 
A time step size At of 2 seconds was'~ chosen initially. 
After 10 sec, At was Increased to 10 sec and to 50 sec after 
200sec. After 1000 sec,At was increased to 100 sec. 


Figure 17 shows the evolution of the temperature of the 
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center interior points of the cylinder and Figure 18 
compares the analytical values with the FEAP’ solution for 
the radial temperature distribuition for two different 


times. 


D. SLAB WITH RADIATION-CONVECTION BOUNDARY CONDITIONS 

Consider a plane slab of thickness L, with constant 
thermophysical properties (k= pe c=1), subjected to 
simultaneous radiation and convection at x=0 and insulated 
at x=L (0<x<L). The slab is initially at a un! form 
temperature Ty and it was decided to consider zero ambient 
temperatures for convection and radiation. 

It was also decided to consider the Biot number 

Bi = h.L/k = 1 
where hes is the convection heat tranfer coefficient and 

R = coT?L/k = 4 
where ¢ is the emissivity and o the Stefan-Boltzman 
constant. 

The same mesh as in Example B was’ employed for’ the 
finite element solution. The time integration was performed 
by the three-time level scheme with y =1/2 and §8=1/3. A 
variable time step size was chosen. A value of At=0.001 was 
used initially. After t=0.02, At was increased to 0.0025. 
After t=0.2, At was increased to 0.01. After t=0.4, At was 
Increased to 0.025. After t=1.0, At was increased to 0.05. 

The history of the ratios T/T) and Tr’ are plotted in 
Figure 19 against the Fourier number 


Fo = at/L- 
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lbs is the surface tmperature at x=0 and (a is the surface 
temperature at x=L. 

The finite element results are compared with those given 
by Haji-Sheick and Sparrow in Ref. 12 and obtained froma 


Monte Carlo method. 
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VI. CONCLUS R 


The computer program in this thesis provides an accurate 
and reliable means for solving a variety of Heat Transfer 
problems. The use of this program and efforts to increase 
its versatility are highly encouraged. 

The capabilities of the program, while quite 
significant, can still be improved. The present version is 
designed for an "in-core" solution technique, which 
restricts the problem size to within the computer core-size. 
The capacity to handle large problems may be increased 


"out of core" technique for solving 


through the use of some 
the system of equations and even for building the mesh data. 
In that case the size of the problems’ treated would be 
restricted only by the availability of external storage 
devices. 

In order to check the mesh input and to easily interpret 
the output, a graphics module must’ be included in the 
program. 

Although the time integration module seems’ to provide 
reliable results in linear and some nonlinear’ problems, it 


should be subjected to further research in order to test it 


with different types of nonlinear analysis. 
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Figure 1. Mesh for Hollow Cylinder Problem 


Figure 2. Eight-Node Serendipity Quadrilateral 
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Figure 3. Local Node Number Sequence for a Quadratic 
Lagrangian Element 
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Figure 4. Local Node Number Sequence for a 21-Node 
Three Dimensional Element 
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Figure 9. (continued) 
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Figure 10. (continued) 
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Figure 11. Step Change in the Loading Term 
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Figure 12. (continued) 
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